Project 7: Mediation Analysis

Causal mediation analysis is frequently used to assess potential causal mechanisms. More specifically, mediation analysis aims at disentangling the effects of a treatment on an outcome through alternative causal mechanisms and has become a popular practice in biomedical and social science applications. The aim of this project is to get familiar with mediation analysis and to understand its main sub-problems: indirect causal effect estimation and sensitivity analysis. You can either choose a more applied approach to apprehend this topic or choose the more theoretical paper on identification, inference and sensitivity.

Our work is based on the following articles:

• LANGE T, VANSTEELANDT S, BEKAERT M. (2011). A Simple Unified Approach for Estimating Natural Direct and Indirect Effects. American Journal of Epidemiology. 176(3):190-195. https://academic.oup.com/aje/article/176/3/190/99496

• IMAI K, KEELE L, YAMAMOTO T. (2010). Identification, Inference and Sensitivity Analysis for Causal Mediation Effects. Statistical Science. 25(1):51-71. https://www.jstor.org/stable/pdf/41058997.pdf?casa_token=ym7l4rAtIZkAAAAA:dDicuRrjqZW-dhh1c_aScoh8OMF7NtX5syaQXE82hYFImMQY_I3UB6DmuD8QpeBpaNRK5fYOxGtgjEad8CVBG4yfVyNzz8m87SYacIUgLOpSDa3HnJkNoQ

• IMAI K, KEELE L, TINGLEY D, YAMAMOTO T. (2010). Causal Mediation Analysis Using R. Advances in Social Science Research Using R. VINOD H. D EDS. Vinod 129-154. https://cran.r-project.org/web/packages/mediation/vignettes/mediation-old.pdf

That Markown is complementray to the Latex report Causal Mediation Analysis and comes after it: we won’t repeat the bibliography for the sake of clarity, as well as the explanations of the methods alrdeady introduced in the report. The Markdown is divided into three parts:

• we first introduce the analysis of the second article based on the package mediation and the dataset JOBS II,
• we then consider a new dataset (European Social Survey) that we pre-process and explore,
• and we eventually apply the methods above on that new dataset.

library(mediation)
library(SparseM)
library(quantreg)
library(nlme)
library(mgcv)
library(foreign)
library(tidyverse)
library(dplyr)
library(rmarkdown)
library(missMDA)
library(stringr)
library(ggplot2)
library(gridExtra)
library(corrplot)
library(cdata)
library(shiny)

1 JOBS II data and mediation package

1.1 Dataset JOBS II

The dataset JOBS II comes from the Job Search Intervention Study (Vinokur and Schul, 1997). In the JOBS II field experiment, 1.801 unemployed workers received a pre-screening questionnaire and were then randomly assigned to treatment and control groups.

Those in the treatment group participated in job-skills workshops while those in the control condition received a booklet describing job-search tips:

• treat: variable descrbing whether the person benefited from the treatment or not.

Two key outcome variables were measured:

• depress2: a continuous measure of post-treatment depressive symptoms based on the Hopkins Symptom Checklist
• work1: a binary variable representing whether the respondent had become employed.

We will successively lead two analysis, (one for each outcome), one variable being continuous and the other discrete (binary).

In the JOBS II data, a key mediating variable is:

• job_seek: a continuous measure of job-search self-efficacy.

Finally, in addition to the outcome and mediators, the JOBS II data also include the following list of baseline covariates that were measured prior to the administration of the treatment:

• depress1: pre-treatment level of depression
• educ: education
• income: income
• nonwhite: race
• marital: marital status
• sex: sex
• occp: previous occupation
• econ_hard: level of economic hardship.

data("jobs")
paged_table(jobs, options = list(rows.print = 10))

1.2 Effects Estimation: Continuous outcome

The study starts with depress2 for the outcome (continuous variable as well as the mediator job_seek).

1.2.1 Linear models: BK procedure

We use the Barron-Kenny procedure to test mediational hypothesis in our causal model: * checking the correlation between the treatment and the outcome, * checking the correlation between the treatment and the mediator, (indeed, no mediation if no effect between the treatment and the mediator), * checking the impact of the correlation between the treatment and the mediator over the outcome.

As explained in our report, the first step is to estimate the regressions coefficients corresponding to the predictions of the mediator M jb_seek and of the outcome Y depress2, with sequential regressions (W referring to the treatment and X to the confounding factors).

• \(Y=cW+\epsilon_1\)

• \(M=aW+\epsilon_2\)

• \(Y={c}'X+bW+\epsilon_3\).

Let’s estimate the coefficients.

model_1_y <- lm(depress2 ~ treat            + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, data = jobs)
model_2_m <- lm(job_seek ~ treat            + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, data = jobs)
model_3_y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, data = jobs)

Here is a summary of the last regression for the outcome.

summary(model_3_y)

## 
## Call:
## lm(formula = depress2 ~ treat + job_seek + depress1 + econ_hard + 
##     sex + age + occp + marital + nonwhite + educ + income, data = jobs)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8247 -0.3834 -0.1090  0.2602  2.8287 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    1.4527281  0.1940639   7.486 1.74e-13 ***
## treat                         -0.0367886  0.0407940  -0.902  0.36740    
## job_seek                      -0.1773802  0.0279535  -6.346 3.56e-10 ***
## depress1                       0.4098612  0.0371003  11.047  < 2e-16 ***
## econ_hard                      0.0679692  0.0221404   3.070  0.00221 ** 
## sex                            0.0621566  0.0448206   1.387  0.16586    
## age                            0.0007858  0.0021942   0.358  0.72034    
## occpmanegerial                 0.0664879  0.0633665   1.049  0.29435    
## occpclerical/kindred           0.0503458  0.0643692   0.782  0.43434    
## occpsales workers             -0.0348333  0.0836727  -0.416  0.67729    
## occpcraftsmen/foremen/kindred -0.0290567  0.0800059  -0.363  0.71656    
## occpoperatives/kindred wrks    0.1635053  0.0829922   1.970  0.04914 *  
## occplaborers/service wrks     -0.0215721  0.0848505  -0.254  0.79937    
## maritalmarried                -0.0072627  0.0551828  -0.132  0.89532    
## maritalsepartd                 0.2019970  0.1132861   1.783  0.07492 .  
## maritaldivrcd                 -0.0453020  0.0639424  -0.708  0.47884    
## maritalwidowed                 0.0923133  0.1458613   0.633  0.52697    
## nonwhitenon.white1            -0.1081444  0.0538549  -2.008  0.04494 *  
## educhighsc                    -0.0023664  0.0901228  -0.026  0.97906    
## educsomcol                     0.0226457  0.0907839   0.249  0.80307    
## educbach                       0.0148269  0.1010784   0.147  0.88341    
## educgradwk                     0.1782504  0.1070053   1.666  0.09611 .  
## income15t24k                  -0.0486597  0.0623912  -0.780  0.43565    
## income25t39k                  -0.0208905  0.0641806  -0.325  0.74488    
## income40t49k                  -0.0528838  0.0780279  -0.678  0.49811    
## income50k+                    -0.1179727  0.0735582  -1.604  0.10912    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.571 on 873 degrees of freedom
## Multiple R-squared:  0.2537, Adjusted R-squared:  0.2324 
## F-statistic: 11.87 on 25 and 873 DF,  p-value: < 2.2e-16

A quick glance at the p-values of the coefficients is enough to ensure that many null hypothesis of non-significance aren’t rejected at \(95\%\): age, marital status… It should also be noted that the the coefficient estimate of pre-level of depression is high (\(0.45\)).

However, we need to be careful in that analysis as the non-significance of a coefficient doesn’t mean that the variable has no impact on the outcome/mediator: that variable could be correlated to another variable concentrating the influence on the predicted variable. Furthermore, correlation does not imply causality AND no correlation does not disaprove causation (Bollen 1989), hence the causal mediation analysis (the causation can be in a mediation relation or hidden by confounding variables).

Here is a summary of the regression for the mediator.

summary(model_2_m)

## 
## Call:
## lm(formula = job_seek ~ treat + depress1 + econ_hard + sex + 
##     age + occp + marital + nonwhite + educ + income, data = jobs)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.84605 -0.39919  0.05924  0.51764  1.46917 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    3.8806256  0.1947174  19.930  < 2e-16 ***
## treat                          0.0774238  0.0492939   1.571 0.116624    
## depress1                      -0.2540256  0.0440638  -5.765 1.13e-08 ***
## econ_hard                      0.1036040  0.0265612   3.901 0.000103 ***
## sex                           -0.0053180  0.0542355  -0.098 0.921913    
## age                            0.0005308  0.0026550   0.200 0.841575    
## occpmanegerial                 0.0056477  0.0766773   0.074 0.941302    
## occpclerical/kindred          -0.1132352  0.0777967  -1.456 0.145883    
## occpsales workers             -0.0137738  0.1012484  -0.136 0.891821    
## occpcraftsmen/foremen/kindred -0.2015647  0.0965720  -2.087 0.037160 *  
## occpoperatives/kindred wrks   -0.2959024  0.0999259  -2.961 0.003147 ** 
## occplaborers/service wrks     -0.3565544  0.1019639  -3.497 0.000494 ***
## maritalmarried                 0.0381422  0.0667623   0.571 0.567934    
## maritalsepartd                 0.3641314  0.1365291   2.667 0.007793 ** 
## maritaldivrcd                  0.2041101  0.0770659   2.649 0.008230 ** 
## maritalwidowed                -0.3301324  0.1761481  -1.874 0.061240 .  
## nonwhitenon.white1             0.0615794  0.0651346   0.945 0.344707    
## educhighsc                     0.1813264  0.1088818   1.665 0.096201 .  
## educsomcol                     0.1638371  0.1097146   1.493 0.135719    
## educbach                       0.2563072  0.1220038   2.101 0.035943 *  
## educgradwk                     0.2013935  0.1293041   1.558 0.119709    
## income15t24k                   0.1583888  0.0753071   2.103 0.035730 *  
## income25t39k                   0.0898314  0.0776033   1.158 0.247355    
## income40t49k                   0.1999402  0.0941763   2.123 0.034031 *  
## income50k+                     0.1631108  0.0888391   1.836 0.066694 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.691 on 874 degrees of freedom
## Multiple R-squared:  0.1243, Adjusted R-squared:  0.1002 
## F-statistic: 5.169 on 24 and 874 DF,  p-value: 2.768e-14

It may seem interesting at first glance to look at the differences with the outcome: the coefficient for pre-level of depression (depress1) is significant and is negative (\(-0.25\)), instead of being positive when predicting the outcome (\(0.4\)), on the other side, variables previously declared as non-significant at \(95\%\) are now declared as significant (marital status). Nonetheless, contrary to basic correlation coefficients, regression coefficients can be affected by the coefficients of the other variables and we had better not rushing to draw conclusions.

Let’s to proceed to the mediation analysis instead. Actually, the package only requires in input the regression models estimated for the second and third equations (model_2_m and model_3_y).

out.1 <- mediate(model_2_m, model_3_y, sims = 1000, treat = "treat", mediator = "job_seek")
summary(out.1)

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME            -0.0140      -0.0317         0.00    0.10
## ADE             -0.0383      -0.1194         0.04    0.34
## Total Effect    -0.0523      -0.1362         0.03    0.21
## Prop. Mediated   0.2210      -1.4844         1.82    0.27
## 
## Sample Size Used: 899 
## 
## 
## Simulations: 1000

The first column correspond to the point estimates and the next two columns to the estimated bounderies of confidence interval at \(95\%\). Here, the confidence intervals are estimated using quasi-Bayesian Monte Carlo approximations (\(1000\) draws). It makes sense as despite the fact that we provide fixed models to the function mediate, the estimation of the effects involve randomness (see equations (11)-(16) in the report).

If we look at the first three columns, it appears that the effects have a very low strength. One can easily confirm that with the the p-values, as they tend to indicate that none of the effects is significant.

Those results are reprensented in the graph below.

plot(out.1)

It appears that those effects (despite being non-significant) are more likely negative; if the job search mediates post-treatment depression, it’s negatively.

Eventually, one can check that - as expected - the Total Effect is the sum of the Average Diret Effect and of the Average Causal Mediation Effect (indirect effect), and that the proportion of the total effect due to mediation corresponds.

1.2.2 Linear models: BK Procedure with Interaction Term

Even if that situation hasn’t been mentioned in the report, it could happen that the treatment moderates the causal mediation effect, i.e. the mediation effect could vary depending on the treatment status. The package mediation enables us to deal with that, using the Baron-Kenny procedure with Interaction term.

Actually, regarding the model it just comes down to consider a kind of mixed effect model, where the coefficients for the prediction of the outcome Y based on the mediator M change with the treatment W. (We can easily do that in R* as we are dealing here with linear models and a binary treatment)*. The model can be written as follows:

model_moder_y <- lm(depress2 ~ treat + job_seek + treat:job_seek +
                        depress1 + econ_hard + sex + age + occp +
                        marital + nonwhite + educ + income, data = jobs)

We display the results below. The intermediate effects have now two potential values, depending on whether the treatment is \(0\) or \(1\). Each time, both values are displayed, in dotted lines for the control group and in solid line for the treatment group.

out.2 <- mediate(model_2_m, model_moder_y, sims = 1000, treat = "treat", mediator = "job_seek")
summary(out.2)

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                          Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control)            -0.0189      -0.0448         0.00    0.10
## ACME (treated)            -0.0120      -0.0281         0.00    0.10
## ADE (control)             -0.0406      -0.1222         0.04    0.32
## ADE (treated)             -0.0337      -0.1130         0.04    0.41
## Total Effect              -0.0526      -0.1359         0.03    0.19
## Prop. Mediated (control)   0.2856      -1.6948         3.29    0.24
## Prop. Mediated (treated)   0.1875      -1.1292         2.24    0.24
## ACME (average)            -0.0154      -0.0351         0.00    0.10
## ADE (average)             -0.0371      -0.1153         0.04    0.35
## Prop. Mediated (average)   0.2365      -1.4329         2.91    0.24
## 
## Sample Size Used: 899 
## 
## 
## Simulations: 1000

plot(out.2, treatment = "both")

The effects don’t seem to be more significant in that framework than in the previous one. It happens that modelling a possible moderation in our situation doesn’t change anything, and it may not be necessary to go into that much details.

1.2.3 Non/Semi-parametric Regression

Considering linear models can sometimes be too restrictive. In that section we investigate options to study mediation effects with non-linear models. For that we consider Generalized Additive Models (GAMs) using splines to model the effects (piecewise polynomial functions).

model_gam_y <- gam(depress2 ~ treat + s(job_seek, bs = "cr") + depress1 + econ_hard +
                     sex + age + occp + marital + nonwhite + educ + income, data = jobs)

Besides, the quasi-Bayesian approximation is suited to linear models for the estimation of the confidence intervals (and more genrally parametric models) , but when it comes to non-linear models bootstrap happens to be the only reliable option (despite being more expensive in terms of computation ressources). We thus use non-parametric bootstrap to estimate confidence intervals in that section, resampling the data with replacement.

out.3 <- mediate(model_2_m, model_gam_y, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")
summary(out.3)

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME            -0.0173      -0.0457         0.00    0.10
## ADE             -0.0264      -0.0995         0.05    0.53
## Total Effect    -0.0437      -0.1212         0.04    0.30
## Prop. Mediated   0.3953      -2.9543         3.78    0.35
## 
## Sample Size Used: 899 
## 
## 
## Simulations: 1000

plot(out.3)

Here no real added-value.

The mediator can also be modelled with a GAM of course, and it is also possible to deal with interaction terms using GAM’s.

1.2.4 Quantile Causal Mediation Effects

Another option that hasn’t been mentioned in the report is to study the distribution of the outcome Y instead of its point values. More specifically, we look at the quantiles of that distribution.

Mediation modeling based on the mean may not be sufficient enough actually to understand different mediating effects across the outcome distribution. For that we replace the previous models for the outcome Y by models of quantile regression, using the package quantreg.

One should notice that the quantile \(tau\) has to be specify in advance. Here we consider the following values for \(tau: 0.2, 0.4, 0.6, 0.8\).

model_quant_y_2 <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, tau = 0.2, data = jobs)

model_quant_y_4 <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, tau = 0.4, data = jobs)

model_quant_y_6 <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, tau = 0.6, data = jobs)

model_quant_y_8 <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, tau = 0.8, data = jobs)

out.4.2 <- mediate(model_2_m, model_quant_y_2, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")
out.4.4 <- mediate(model_2_m, model_quant_y_4, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")
out.4.6 <- mediate(model_2_m, model_quant_y_6, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")
out.4.8 <- mediate(model_2_m, model_quant_y_8, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")

Here is the summary of the effects for \(tau=0.2\).

summary(out.4.2)

## 
## Causal Mediation Analysis 
## 
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
## 
##                Estimate 95% CI Lower 95% CI Upper p-value  
## ACME            -0.0121      -0.0283         0.00   0.084 .
## ADE              0.0419      -0.0494         0.08   0.550  
## Total Effect     0.0298      -0.0590         0.07   0.838  
## Prop. Mediated  -0.4079      -6.5514         6.66   0.890  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 899 
## 
## 
## Simulations: 1000

The plots are grouped below.

par(mfrow=c(2,2))
plot(out.4.2, main="Tau = 0.2")
plot(out.4.4, main="Tau = 0.4")
plot(out.4.6, main="Tau = 0.6")
plot(out.4.8, main="Tau = 0.8")

Only the quantile \(tau = 0.2\) differs from the others, the other results are quite the same as for the estimation with the mean. More specifically, the ACME doesn’t change and the only difference is with the ADE. Even if the ADE effect don“t appear to be significant and the x-axis of the plot \(tau = 0.2\) is more tightened than for the other plots (be careful with that), that effect is more positive than negative with \(tau = 0.2\).

One should notice that this result could hardly be attributed to a possible high variance of the effects for the tails of the distribution, provided that:

• the ACME isn’t affected at all, only the ADE

• the distribution on the left is quite flat, and the quantile \(tau=0.2\) may not be that instable

• the other tail (way thiner) doesn’t face that issue.

The figure below may help support those arguments.

ggplot(jobs) + aes(depress2) + geom_density(alpha=0.6, color="darkolivegreen1", fill="darkolivegreen1") +
   ggtitle("Symptoms of depression, post-treatment", subtitle="Density distribution")

We can thus affirm that - even if the effect isn’t significant - the treatment is more likely to increase depressive symptoms for the people almost not depressed at all than for the others. Nonetheless, one must keep in mind that those effects aren’t significant and that the confidence intervals don’t enable us to go further.

Apart from that quantile, we can nonetheless say that the mediating effects don’t vary much across the distribution of the outcome.

NB: we could also look at the distribution of the difference depress2 - depress1, symptoms of depression pre and post treatment (see below), as depress1 is taken into account in the regression.

ggplot(jobs) + aes(x=jobs$depress2 - jobs$depress1) +
  geom_density(alpha=0.6, color="darkolivegreen1", fill="darkolivegreen1") +
  scale_x_discrete(name="depress2 - depress1", limits=c(0:6)-3) +
  ggtitle("Difference of symptoms of depression, pre and post treatment", subtitle="Density distribution")

1.3 Effects Estimation: Discrete outcome

Discrete outcomes can’t be treated as continuous outcomes (for the regressions amongs other things). Nonetheless, the package mediation handles that issue. We illustrate that with the discrete outcome work1.

To find out the estimate we use a probit regression instead of a linear regression. We could choose a logit regression, but the sensitivity analysis of the package mediation requires a probit regression.

model_discrete_y <- glm(work1 ~ treat + job_seek + depress1 + econ_hard + sex + age +
                                occp + marital + nonwhite + educ + income,
                                family = binomial(link = "probit"), data = jobs)

out.5 <- mediate(model_2_m, model_discrete_y, sims = 1000, treat = "treat", mediator = "job_seek")
summary(out.5)

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                           Estimate 95% CI Lower 95% CI Upper p-value  
## ACME (control)            0.003149    -0.000933         0.01   0.146  
## ACME (treated)            0.003387    -0.000982         0.01   0.146  
## ADE (control)             0.055884    -0.001460         0.11   0.064 .
## ADE (treated)             0.056123    -0.001468         0.11   0.064 .
## Total Effect              0.059272     0.000880         0.11   0.050 *
## Prop. Mediated (control)  0.045554    -0.052445         0.44   0.184  
## Prop. Mediated (treated)  0.049405    -0.053708         0.45   0.184  
## ACME (average)            0.003268    -0.000964         0.01   0.146  
## ADE (average)             0.056004    -0.001464         0.11   0.064 .
## Prop. Mediated (average)  0.047480    -0.052494         0.44   0.184  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 899 
## 
## 
## Simulations: 1000

plot(out.5)

Although we have no interaction term in our model, we still have to a difference between the treatment and control groups, the gap can be overlooked as simply due to the non-linearity of the model.

1.4 Sensitivity Analysis

Here , the idea of the sensitivity analysis is to study the robustness of our results to the ignorability assumption (and more precisely to the second condition as explained in the report). Indeed, conclusions that are harder to alter by a plausible violation of an identification assumption add a higher scientific value to analytics.

We keep considering \(1000\) simulations per call to the function mediate - as before - due to the variance of our results that require a high number of simulations (we have tried samller values but that wasn’t convincing enough).

Once again, we use the Baron-Kenny Procedure.

model_2_m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age + occp +
                  marital + nonwhite + educ + income, data = jobs)

model_3_y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp +
                  marital + nonwhite + educ + income, data = jobs)

As in the report, we present two ways of leading the sensitivity analysis: through the study of the residuals correlation \(\rho\), or through the coefficients of determination \(R\). The mediation package can do both.

med_cont <- mediate(model_2_m, model_3_y, sims=1000, treat = "treat", mediator = "job_seek")

sens_cont <- medsens(med_cont, rho.by = 0.05)

The rho.by option neables us to choose the precision for the incrementation of the parameter rho in the sensitivity analysis. Having a rough incrementation goes faster but don’t allow us to estimate precisely the robustness of our results. We choose the same value of \(\rho\) as in the article.

The object sens_cont contains the estimated effects and their CI’s corresponding to each value of \(\rho\) and \(R\).

summary(sens_cont)

## 
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
## 
## Sensitivity Region
## 
##         Rho    ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
##  [1,] -0.95  0.1808      -0.0418       0.4034       0.9025       0.5898
##  [2,] -0.90  0.1183      -0.0274       0.2640       0.8100       0.5293
##  [3,] -0.85  0.0894      -0.0207       0.1996       0.7225       0.4722
##  [4,] -0.80  0.0715      -0.0166       0.1597       0.6400       0.4182
##  [5,] -0.75  0.0588      -0.0137       0.1312       0.5625       0.3676
##  [6,] -0.70  0.0489      -0.0115       0.1093       0.4900       0.3202
##  [7,] -0.65  0.0410      -0.0096       0.0916       0.4225       0.2761
##  [8,] -0.60  0.0342      -0.0081       0.0766       0.3600       0.2353
##  [9,] -0.55  0.0284      -0.0068       0.0636       0.3025       0.1977
## [10,] -0.50  0.0232      -0.0057       0.0520       0.2500       0.1634
## [11,] -0.45  0.0185      -0.0047       0.0416       0.2025       0.1323
## [12,] -0.40  0.0142      -0.0038       0.0321       0.1600       0.1046
## [13,] -0.35  0.0102      -0.0030       0.0233       0.1225       0.0801
## [14,] -0.30  0.0064      -0.0025       0.0153       0.0900       0.0588
## [15,] -0.25  0.0028      -0.0026       0.0082       0.0625       0.0408
## [16,] -0.20 -0.0007      -0.0049       0.0036       0.0400       0.0261
## [17,] -0.15 -0.0040      -0.0105       0.0025       0.0225       0.0147
## [18,] -0.10 -0.0073      -0.0172       0.0026       0.0100       0.0065
## [19,] -0.05 -0.0105      -0.0242       0.0031       0.0025       0.0016
## [20,]  0.00 -0.0137      -0.0312       0.0037       0.0000       0.0000
## [21,]  0.05 -0.0169      -0.0382       0.0043       0.0025       0.0016
## [22,]  0.10 -0.0202      -0.0453       0.0050       0.0100       0.0065
## [23,]  0.15 -0.0234      -0.0526       0.0057       0.0225       0.0147
## [24,]  0.20 -0.0268      -0.0600       0.0065       0.0400       0.0261
## [25,]  0.25 -0.0302      -0.0677       0.0072       0.0625       0.0408
## [26,]  0.30 -0.0338      -0.0757       0.0080       0.0900       0.0588
## [27,]  0.35 -0.0376      -0.0841       0.0089       0.1225       0.0801
## [28,]  0.40 -0.0416      -0.0931       0.0098       0.1600       0.1046
## [29,]  0.45 -0.0460      -0.1027       0.0108       0.2025       0.1323
## [30,]  0.50 -0.0507      -0.1132       0.0118       0.2500       0.1634
## [31,]  0.55 -0.0558      -0.1247       0.0130       0.3025       0.1977
## [32,]  0.60 -0.0617      -0.1378       0.0144       0.3600       0.2353
## [33,]  0.65 -0.0684      -0.1528       0.0159       0.4225       0.2761
## [34,]  0.70 -0.0764      -0.1706       0.0177       0.4900       0.3202
## [35,]  0.75 -0.0862      -0.1925       0.0200       0.5625       0.3676
## [36,]  0.80 -0.0990      -0.2209       0.0229       0.6400       0.4182
## [37,]  0.85 -0.1169      -0.2609       0.0271       0.7225       0.4722
## [38,]  0.90 -0.1458      -0.3252       0.0337       0.8100       0.5293
## [39,]  0.95 -0.2083      -0.4647       0.0482       0.9025       0.5898
## 
## Rho at which ACME = 0: -0.2
## R^2_M*R^2_Y* at which ACME = 0: 0.04
## R^2_M~R^2_Y~ at which ACME = 0: 0.0261

1.4.1 Residuals correlation

plot(sens_cont, sens.par = "rho")

We recall that the model is identifiable supposing that we know the correlation \(\rho\) between \(\epsilon_2\) and \(\epsilon_3\) (see report and part 1.2.1 of the Markdown).

Let’s explain what is depicted on the graph.

• The solid line curve represents the mediation effect ACME against different values of \(\rho\), with the \(95%\) confidence interval (based on the delta method).

• The dashed horizontal line represents the estimated mediation effect under the sequential ignorability assumption (previous quantification of the effect). It thus correspond to the effect we get when considering \(\rho = 0\) (with our current model).

• The two extremities of the graph correspond to \(\rho = +/- 1\), i.e. a total correlation of \(\epsilon_2\) and \(\epsilon_3\) and thus (see part 1.2.1 and report) a total correlation of \(Y_i(w', m)\) and \(M(w)\) conditional to \(W=w\) (for any \(w, w', m\)), hinting thus at a very high mediation effect.

• Why an \(AMCE\) positive when \(\rho\) becomes low enough (closer to \(-1\)) ? Because it corresponds to an unobserved pre-treatment covariate(s) that makes \(Y\) decrease when it makes \(M\) increase. By adjusting with regards to that covariate we would get that \(Y\) increases more (or decreases less) when \(M\) increases (through the intermediate effect).

How to interpret it ?

The critical point is when the curve crosses the x-axis: the mediation effect nullifies. That point corresponds here to \(\rho = -0.20\) (see summary above). Moreover, we can see that the \(95\%\) confidence intervals always contain \(0\), whatever the value of \(\rho\), and that our conclusion of non-significance of the \(AMCE\) effect may hold with any value of \(\rho\).

If we had had a significant effect previously, we would have looked at the range of \(\rho\) on which the \(95\%\) confidence interval \(CI\) don’t contain \(0\) in order to determine whether the results are robust. If that condition of \(0 \notin CI\) would only be true for a small interval of \(\rho\) around \(0\), that would be disappointed, otherwise it would hint at the fact that our results resist to a violation of the second condition required for ignorability.

1.4.2 Coefficients of determination

We look at the products of the coefficients \(R^{2*}_M R^{2*}_Y\), corresponding to the proportion of the variance of our outcome which is unexplained by our current predictors, and that would be explained by a non-identified pre-treatment covariate (see report).

By the way, we have \(\rho^2 = R_M^{2*}R_Y^{2*}\) (see report) and \(\rho = sign(\lambda_2 \lambda_3) R_M^* R_Y^*\) with \(sign(\lambda_2 \lambda_3)\) the sign of the product; positive if the correlation with the unobserved covariate is of the same sign for the outcome and the mediator (see report), negative otherwise.

As above, the idea is to look at the estimated value of the \(AMCE\) for different values of \(R^{2*}_M\) and \(R^{2*}_Y\) that we choose, “what if the unobserved covariate had that precise importance ?”. More specifically the plot represents level curves of the \(AMCE\), (when the product is set).

par(mfrow=c(1,2))
plot(sens_cont, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens_cont, sens.par = "R2", r.type = "residual", sign.prod = "positive")

The level curve corresponding to \(ACME = 0\) is in bold (critical point in the plot with \(\rho\)).

It is important to understand that the shape of the curves will always be the same as they correspond to \(xy= c^{te}\), for different constant values. (There are two plots to study both signs, as the two plots may not give the same values of ACME).

The point is thus to look at the value of ACME on each level curve: does our result hold when the product increases ? i.e. when we consider a neglected pre-treatment covariate with more influence, do we still get that the effect is significant ? The signs for extrem values can be explained as with \(\rho\).

This is only an graphical visualization of the summary, helping us to visualize when the ACME wouldn’t be really significant anymore.

Here the effect was already non-significant, and even with different values of the product it seems to remain non-significant. Nonetheless, if we originally had a significant effect, let’s say negative, we would have looked for the value of the product where the effect isn’t strictly negative anymore. The summary above can help. That would have indicated us the proportion of variance explained by neglected confounders which is required to have a situation where the effect would actually be null. That would help us having more hindsight on the robustness of our results by determining whether they would resist a violation of the ignorability assumption.

2 European Social Survey

2.1 Presenting the dataset

The new data we chose come from the European Social Survey (ESS), a European initiative to study the attitudes, beliefs and behaviour patterns of various populations in Europe. The data come from a series of questionnaires collected every 2 years since 2002. They are available for free at the following link: https://www.europeansocialsurvey.org/data#.

The list of all the available variables is given here : https://www.europeansocialsurvey.org/docs/cumulative/ESS_cumulative_variable_list.pdf. One may have noticed that the the questionnaires change from one version to another, the study covering more and more topics over time, and some variables are thus only available for some precise years.

Our work focuses on studying opinions on immigration questions: we thus had to choose between the first round (2002) and the seventh round (2014) of the survey. We chose the most recent round, available here: https://www.europeansocialsurvey.org/data/download.html?r=7. Those version contain indeed questions measuring attitudes of immigration and perceptions of social realities as well as opinions on public policy and knowledge, in addition to the standard questions that don’t vary over the questionnaires.

The table below contains the questions we decided to select from the raw data (40185 rows).

d <- read.dta("D:/Téléchargements/ESS7e02_2.stata/ESS7e02_2.dta")

# columns we keep
df <- d %>% dplyr::select(c(acetalv, imbleco, imtcjob, imwbcrm, imueclt, gvrfgap, imwbcnt, gndr, yrbrn, eduyrs, uempla, stfeco, tvtot, tvpol, happy, sclmeet, stflife, rlgblg, dscrgrp, ctzcntr, brncntr, blgetmg, facntr, mocntr, ipeqopt, impsafe, iphlppl, ipstrgv, pplstrd))

paged_table(df, options = list(rows.print = 10))

2.2 Presenting our study

Our objective is to decompose the thinking process leading to a specific opinion on immigration by breaking it down into intermediate steps.

We try to better explain basic political opinions on immigration by the fact that minority race/ethnic groups live in the same local area or not. Our work aims specifically at determining whether those opinions involve intermediate (mediating) ideas such as “immigrants would take jobs away”, “worsen criminality”, “undermine the country’s cultural life”…, as well as the influence of living close to minorities on those ideas.

Are such ideas more present in the minds of people living close to minorities, or on the contrary in people having no real contact with them ? potentially hinting at the fact that some ideas could actually be prejudices, or not.

One thing to note is that we don’t study all those mediating ideas (effects) at the same time, only one after another (see the report fore more details regarding the issue of the multiplicity of mediators). Same for the outcomes.

We also adjust by gender, age, education, job situation, time watching tv, belonging to a minority and many other confounding factors listed below.

Here are the questions we consider in our study:

EXPOSURE W

acetalv: People of minority race/ethnic group in current living area

MEDIATORS M (choose one)

imbleco: Taxes and services: immigrants take out more than they put in or less

imtcjob: Immigrants take jobs away in country or create new jobs

imwbcrm: Immigrants make country’s crime problems worse or better

imueclt: Country’s cultural life undermined or enriched by immigrants

OUTCOME Y (choose one)

gvrfgap: Government should be generous judging applications for refugee status

imwbcnt: Immigrants make country worse or better place to live

CONFOUNDING FACTORS X

gndr: Gender

yrbrn: Year of birth

eduyrs: Years of full-time education completed

uempla: Doing last 7 days: unemployed, actively looking for job

stfeco: How satisfied with present state of economy in country

tvtot: TV watching, total time on average weekday

tvpol: TV watching, news/ politics/current affairs on average weekday

happy: How happy are you

sclmeet: How often socially meet with friends, relatives or colleagues

stflife: How satisfied with life as a whole

rlgblg: Belonging to particular religion or denomination

dscrgrp: Member of a group discriminated against in this country

ctzcntr: Citizen of country

brncntr: Born in country

blgetmg: Belong to minority ethnic group in country

facntr: Father born in country

mocntr: Mother born in country

ipeqopt: Important that people are treated equally and have equal opportunities

impsafe: Important to live in secure and safe surroundings

iphlppl: Important to help people and care for others well-being

ipstrgv: Important that government is strong and ensures safety

pplstrd: Better for a country if almost everyone share customs and traditions.

2.3 Pre-processing

2.3.1 Categories values

We attribute a numeric value to every answer (instead of a textual categorical value), based on the information provided by the Norwegian Social Science Data Services (NSD). Most of the variables are indeed ordinal variables but not in a numeric format yet.

# Replacement functions

rep_no_ans <- function(x) {
  if (is.na(x)) NaN
  else if (x %in% c("Refusal", "Don't know", "No answer", "Not applicable")) NaN
  else x}

rep_acetalv <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Almost nobody minority race/ethnic group") 0
  else if (x == "Some minority race/ethnic group") 1
  else if (x == "Many minority race/ethnic group") 2
  else NaN}

rep_imbleco <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Generally take out more") 0
  else if (x == "Generally put in more") 10
  else x}

rep_imtcjob <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Take jobs away") 0
  else if (x == "Create new jobs") 10
  else x}

rep_imwbcrm <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Crime problems made worse") 0
  else if (x == "Crime problems made better") 10
  else x}

rep_imueclt_happy <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Cultural life undermined") 0
  else if (x == "Cultural life enriched") 10
  else x}

rep_agree <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Agree strongly") 1
  else if (x == "Agree") 2
  else if (x == "Neither agree nor disagree") 3
  else if (x == "Disagree") 4
  else if (x == "Disagree strongly") 5
  else x}

rep_imwbcnt <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Worse place to live") 0
  else if (x == "Better place to live") 10
  else x}

rep_gndr <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Male") 0
  else if (x == "Female") 1
  else x}

rep_uempla <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Not marked") 0
  else if (x == "Marked") 1
  else x}

rep_stf <- function(x) {
  if (is.na(x)) NaN
  else if (x == "   Extremely dissatisfied") 0
  else if (x == "   Extremely satisfied") 10
  else x}

rep_tv <- function(x) {
  if (is.na(x)) NaN
  else if (x == "No time at all") 0
  else if (x == "Less than 0,5 hour") 1
  else if (x == "0,5 hour to 1 hour") 2
  else if (x == "More than 1 hour, up to 1,5 hours") 3
  else if (x == "More than 1,5 hours, up to 2 hours") 4
  else if (x == "More than 2 hours, up to 2,5 hours") 5
  else if (x == "More than 2,5 hours, up to 3 hours") 6
  else if (x == "More than 3 hours") 7
  else x}

rep_sclmeet <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Never") 1
  else if (x == "Less than once a month") 2
  else if (x == "Once a month") 3
  else if (x == "Several times a month") 4
  else if (x == "Once a week") 5
  else if (x == "Several times a week") 6
  else if (x == "Every day") 7
  else x}

rep_yes_no <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Yes") 0
  else if (x == "No") 1
  else x}

rep_like <- function(x) {
  if (is.na(x)) NaN
  else if (x == "Very much like me") 1
  else if (x == "Like me") 2
  else if (x == "Somewhat like me") 3
  else if (x == "A little like me") 4
  else if (x == "Not like me") 5
  else if (x == "Not like me at all") 6
  else x}

# Applying those functions to the dataset

df$acetalv <- Vectorize(rep_acetalv)(df$acetalv)
df$imbleco <- Vectorize(rep_imbleco)(df$imbleco)
df$imtcjob <- Vectorize(rep_imtcjob)(df$imtcjob)
df$imwbcrm <- Vectorize(rep_imwbcrm)(df$imwbcrm)
df$imueclt <- Vectorize(rep_imueclt_happy)(df$imueclt)
df$gvrfgap <- Vectorize(rep_agree)(df$gvrfgap)
df$imwbcnt <- Vectorize(rep_imwbcnt)(df$imwbcnt)
df$gndr    <- Vectorize(rep_gndr)(df$gndr)
df$uempla  <- Vectorize(rep_uempla)(df$uempla)
df$stfeco  <- Vectorize(rep_stf)(df$stfeco)
df$tvtot   <- Vectorize(rep_tv)(df$tvtot)
df$tvpol   <- Vectorize(rep_tv)(df$tvpol)
df$happy   <- Vectorize(rep_imueclt_happy)(df$happy)
df$sclmeet <- Vectorize(rep_sclmeet)(df$sclmeet)
df$stflife <- Vectorize(rep_stf)(df$stflife)
df$rlgblg  <- Vectorize(rep_yes_no)(df$rlgblg)
df$dscrgrp <- Vectorize(rep_yes_no)(df$dscrgrp)
df$ctzcntr <- Vectorize(rep_yes_no)(df$ctzcntr)
df$brncntr <- Vectorize(rep_yes_no)(df$brncntr)
df$facntr  <- Vectorize(rep_yes_no)(df$facntr)
df$mocntr  <- Vectorize(rep_yes_no)(df$mocntr)
df$ipeqopt <- Vectorize(rep_like)(df$ipeqopt)
df$impsafe <- Vectorize(rep_like)(df$impsafe)
df$iphlppl <- Vectorize(rep_like)(df$iphlppl)
df$ipstrgv <- Vectorize(rep_like)(df$ipstrgv)
df$pplstrd <- Vectorize(rep_agree)(df$pplstrd)


df <- df %>% dplyr::mutate_all(Vectorize(rep_no_ans))
df <- df %>% dplyr::mutate_all(as.numeric)
df <- df %>% dplyr::mutate_at(c("gndr", "uempla", "rlgblg", "dscrgrp", "ctzcntr", "brncntr", "blgetmg", "facntr", "mocntr"), as.factor)
df <- replace(df, df =="NaN", NA)

2.3.2 Imputing missing data

We impute the missing data with an iterative Factorial Analysis for Mixed Data model (FAMD) imputation method.

Audigier, V., Husson, F. & Josse, J. (2013). A principal components method to impute mixed data. Advances in Data Analysis and Classification, 10(1), 5-26. https://arxiv.org/abs/1301.4797

df_imputed <- imputeFAMD(df, ncp=10)$completeObs
df_imputed <- df_imputed %>% mutate_if(is.factor, function(x) as.factor(str_sub(x,-1))) %>% mutate_if(is.numeric, round)

The final (clean) dataset:

paged_table(df_imputed, options = list(rows.print = 10))

2.4 Data exploration

2.4.1 Outcomes

As a recall, the outcome variables are:

• gvrfgap: Government should be generous judging applications for refugee status (1: Agree strongly - 5: Disagree strongly)

• imwbcnt: Immigrants make country worse or better place to live (0: Worse - 10: Better).

Here are histograms for both variables, with the imputed data in red and the raw data in blue.

ggplot() + geom_bar(aes(x=df_imputed$gvrfgap), alpha=0.2, color="red",  fill="red") +
           geom_bar(aes(x=df$gvrfgap        ), alpha=0.2, color="blue", fill="blue") +
           scale_x_discrete(name ="gvrfgap", limits=c(1:5)) +
           labs(title = "Government should be generous judging applications for refugee status",
                subtitle = "(1: Agree strongly - 5: Disagree strongly)",
                caption = "Data source: ESS Round 7")

ggplot() + geom_bar(aes(x=df_imputed$imwbcnt), alpha=0.2, color="red",  fill="red") +
           geom_bar(aes(x=df$imwbcnt        ), alpha=0.2, color="blue", fill="blue") +
           scale_x_discrete(name ="imwbcnt", limits=c(0:10)) +
           labs(title = "Immigrants make country worse or better place to live",
                subtitle = "(0: Worse - 10: Better)",
                caption = "Data source: ESS Round 7")

Sounds like we have values Missing not at Random (MNAR) for the outcomes: with no surprise it seems - based on the imputation results - that the people not answering are more likely to have a neutral opinion (are at least look like the people having answered with a neutral opinion).

We expect the two potential outcomes to be negatively correlated, let’s compare them.

ggplot(df_imputed) + aes(x = gvrfgap, y=imwbcnt) + geom_jitter(color="orange", size=0.07) +
             scale_y_discrete(limits=c(0:10)) +
             labs(title = "Graphical correlation of the potential outcomes",
                  subtitle = "Jitter plot",
                  caption = "Data source: ESS Round 7")

We choose the outcome question gvrfgap: “Government should be generous judging applications for refugee status”.

2.4.2 Exposure

The exposure is acetalv: People of minority race/ethnic group in current living area.

Imputed data in orange and non-imputed data in green.

ggplot() + geom_bar(aes(x=df_imputed$acetalv), alpha=0.2, color="orange", fill="orange") +
           geom_bar(aes(x=df$acetalv        ), alpha=0.2, color="green", fill="green") +
           scale_x_discrete(name ="acetalv", limits=c(0:2)) +
           labs(title = "People of minority race/ethnic group in current living area",
                subtitle = "(0: Almost nobody - 2: Many minority groups)",
                caption = "Data source: ESS Round 7")

2.4.3 Mediators

The list of possible mediators is:

• imbleco: Taxes and services: immigrants take out more than they put in or less

• imtcjob: Immigrants take jobs away in country or create new jobs

• imwbcrm: Immigrants make country’s crime problems worse or better

• imueclt: Country’s cultural life undermined or enriched by immigrants.

mediators = df_imputed %>% select(imbleco, imtcjob, imwbcrm, imueclt)

mediators <- sample_n(mediators, 2000)

meas_vars <- c("imbleco", "imtcjob", "imwbcrm", "imueclt")

controlTable <- data.frame(expand.grid(meas_vars, meas_vars, stringsAsFactors = FALSE))
colnames(controlTable) <- c("xv", "yv")
controlTable <- cbind(data.frame(pair_key = paste(controlTable[[1]], controlTable[[2]]),
                      stringsAsFactors = FALSE), controlTable)
aug = rowrecs_to_blocks(mediators, controlTable)
splt <- strsplit(aug$pair_key, split = " ", fixed = TRUE)
aug$x <- vapply(splt, function(si) si[[1]], character(1))
aug$y <- vapply(splt, function(si) si[[2]], character(1))
aug$x <- factor(as.character(aug$x), meas_vars)
aug$y <- factor(as.character(aug$y), meas_vars)


ggplot(aug, aes(x=xv, y=yv)) +
       geom_jitter(color="darkslategray1", size=0.001) +
       scale_color_brewer("Diamond\nclarity") +
       facet_grid(y~x, labeller = label_both, scale = "free") +
       ggtitle("Mediators") +
       ylab(NULL) + xlab(NULL) +
       theme_dark()

As we consider small plots (having many plots in the same graph), a jitter plot can’t be refined enough to distinguish between dense areas, hence the density plot in 2D below. One may note that the areas with few points appear as empty.

The mediators tend to be pretty correlated, we analyze them one at at time.

2.4.4 Confounding factors (Shiny app)

We offer an app to visualize any confounder and its relation with any “main variable”" (outcome, exposure, mediator).

On the left plot (confounder only) the blue light portion conrrespond to the missing data we imputed.

It is possible to choose the type of graph to study correlations between confounder / main variable.

## PhantomJS not found. You can install it with webshot::install_phantomjs(). If it is installed, please make sure the phantomjs executable can be found via the PATH variable.

Shiny applications not supported in static R Markdown documents

3 Mediation Analysis on the new dataset

The objective of that part is to conduct a mediation analysis on our new dataset. We will use the methods explain in part 1.

We will consider the treatment acetalv, but as binary treatment: the groups “Some minorities” and “Many minorities” will thus be grouped together.

dd <- df_imputed %>% mutate_all(as.numeric)
dd$acetalv[dd$acetalv==2] <- 1

The outcome we chose is gvrfgap, but we will test all the mediators: imbleco, imtcjob, imwbcrm and imueclt. We select some of the confounding factors introduced below.

We quantify the effect with linear models, as non-linear models don’t bring any change (not shown here but we’ve looked at it). We have also looked at interaction terms and it makes very little difference too.

3.1 Estimations and tests of the effects

# mediator imbleco

model.m.1 <- lm(imbleco ~ acetalv           + gndr + yrbrn + uempla + stfeco + tvtot +
                  happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)

model.y.1 <- lm(gvrfgap ~ acetalv + imbleco + gndr + yrbrn + uempla + stfeco + tvtot +
                  happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)

out.1 <- mediate(model.m.1, model.y.1, sims = 1000, treat = "acetalv", mediator = "imbleco")

# mediator imtcjob

model.m.2 <- lm(imtcjob ~ acetalv           + gndr + yrbrn + uempla + stfeco + tvtot +
                  happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)

model.y.2 <- lm(gvrfgap ~ acetalv + imtcjob + gndr + yrbrn + uempla + stfeco + tvtot +
                  happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)

out.2 <- mediate(model.m.2, model.y.2, sims = 1000, treat = "acetalv", mediator = "imtcjob")

# mediator imwbcrm

model.m.3 <- lm(imwbcrm ~ acetalv           + gndr + yrbrn + uempla + stfeco + tvtot +
                  happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)

model.y.3 <- lm(gvrfgap ~ acetalv + imwbcrm + gndr + yrbrn + uempla + stfeco + tvtot +
                  happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)

out.3 <- mediate(model.m.3, model.y.3, sims = 1000, treat = "acetalv", mediator = "imwbcrm")

# mediator imueclt

model.m.4 <- lm(imueclt ~ acetalv           + gndr + yrbrn + uempla + stfeco + tvtot +
                  happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)

model.y.4 <- lm(gvrfgap ~ acetalv + imueclt + gndr + yrbrn + uempla + stfeco + tvtot +
                  happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)

out.4 <- mediate(model.m.4, model.y.4, sims = 1000, treat = "acetalv", mediator = "imueclt")

We plot the results of the four mediators together.

par(mfrow=c(2,2))
plot(out.1, main="imbleco")
plot(out.2, main="imtcjob")
plot(out.3, main="imwbcrm")
plot(out.4, main="imueclt")

When looking at the summaries below, we can see that the three effects are almost always significant at \(95\%\) (one can look at the p-values but also at the confidence intervals on the plots that don’t contain \(0\)).

The Total Effect is relatively stable and the proportion of mediated effect really depends on the mediator.

We have for:

• imbleco around \(10\%\)
• imtcjob around \(25-30\%\)
• imwbcrm around \(5\%\)
• imueclt around \(30\%\).

The rest is explained by the Direct Effect.

Based on the proportion of mediated effect in the total effect, we can distinguish between two groups of mediators:

• imbleco and imwbcrm on the left: low mediation

• imtcjob and imueclt on the right: explain almost one third of the total effect.

Here is the summary for imbleco: Taxes and services: immigrants take out more than they put in or less.

summary(out.1)

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME           -0.00935     -0.01496         0.00  <2e-16 ***
## ADE            -0.08498     -0.10511        -0.06  <2e-16 ***
## Total Effect   -0.09432     -0.11577        -0.07  <2e-16 ***
## Prop. Mediated  0.09827      0.04271         0.16  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 40185 
## 
## 
## Simulations: 1000

Here is the summary for imtcjob: Immigrants take jobs away in country or create new jobs.

summary(out.2)

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME            -0.0259      -0.0311        -0.02  <2e-16 ***
## ADE             -0.0680      -0.0877        -0.05  <2e-16 ***
## Total Effect    -0.0939      -0.1147        -0.07  <2e-16 ***
## Prop. Mediated   0.2770       0.2140         0.37  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 40185 
## 
## 
## Simulations: 1000

Here is the summary for imwbcrm: Immigrants make country’s crime problems worse or better.

summary(out.3)

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME           -0.00443     -0.00953         0.00   0.076 .  
## ADE            -0.09002     -0.11038        -0.07  <2e-16 ***
## Total Effect   -0.09445     -0.11472        -0.07  <2e-16 ***
## Prop. Mediated  0.04610     -0.00477         0.10   0.076 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 40185 
## 
## 
## Simulations: 1000

Here is the summary for imueclt: Country’s cultural life undermined or enriched by immigrants.

summary(out.4)

## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                Estimate 95% CI Lower 95% CI Upper p-value    
## ACME            -0.0293      -0.0357        -0.02  <2e-16 ***
## ADE             -0.0657      -0.0847        -0.05  <2e-16 ***
## Total Effect    -0.0949      -0.1162        -0.07  <2e-16 ***
## Prop. Mediated   0.3094       0.2373         0.39  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 40185 
## 
## 
## Simulations: 1000

3.2 Sensitivity analysis

We conduct the same kind of sensitivity analysis as with the dataset JOBS II in part 1, i.e looking at the residual coefficient \(\rho\) and at the coefficients of determination \(R^*\).

med.cont.1 <- mediate(model.m.1, model.y.1, sims=1000, treat = "acetalv", mediator = "imbleco")
med.cont.2 <- mediate(model.m.2, model.y.2, sims=1000, treat = "acetalv", mediator = "imtcjob")
med.cont.3 <- mediate(model.m.3, model.y.3, sims=1000, treat = "acetalv", mediator = "imwbcrm")
med.cont.4 <- mediate(model.m.4, model.y.4, sims=1000, treat = "acetalv", mediator = "imueclt")

sens.cont.1 <- medsens(med.cont.1, rho.by = 0.05)
sens.cont.2 <- medsens(med.cont.2, rho.by = 0.05)
sens.cont.3 <- medsens(med.cont.3, rho.by = 0.05)
sens.cont.4 <- medsens(med.cont.4, rho.by = 0.05)

3.2.1 Residual coefficient

The residuals correlation plots are grouped below:

par(mfrow=c(2,2))
plot(sens.cont.1, sens.par = "rho", sub="imbleco")
plot(sens.cont.2, sens.par = "rho", sub="imtcjob")
plot(sens.cont.3, sens.par = "rho", sub="imwbcrm")
plot(sens.cont.4, sens.par = "rho", sub="imueclt")

With no surprise we can see that the higher the proportion of mediated effect, the more robust the result: indeed, for the mediators on the right, the effect remains strictly negative at \(95\%\) for a pretty wide interval of values \(\rho\) (the confidence intervals are also very thin). The confidence intervals for the first mediator imbleco are less thin but the effect remains strictly negative at \(95\%\) on a wide interval. The real issue regarding robustness concerns the mediator imwbcrm: the confidence intervals are large and contain zero for a a wide range of \(\rho\) values.

Besides, we can see on the summaries (below) that the ACME nullifies at \(\rho \leq -0.25\) for all mediators.

Here is the summary for imbleco: Taxes and services: immigrants take out more than they put in or less.

summary(sens.cont.1)

## 
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
## 
## Rho at which ACME = 0: -0.25
## R^2_M*R^2_Y* at which ACME = 0: 0.0625
## R^2_M~R^2_Y~ at which ACME = 0: 0.0494

Here is the summary for imtcjob: Immigrants take jobs away in country or create new jobs.

summary(sens.cont.2)

## 
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
## 
## Sensitivity Region
## 
##        Rho  ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
## [1,] -0.25 9e-04       -1e-04       0.0019       0.0625       0.0479
## 
## Rho at which ACME = 0: -0.25
## R^2_M*R^2_Y* at which ACME = 0: 0.0625
## R^2_M~R^2_Y~ at which ACME = 0: 0.0479

Here is the summary for imwbcrm: Immigrants make country’s crime problems worse or better.

summary(sens.cont.3)

## 
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
## 
## Sensitivity Region
## 
##         Rho    ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
##  [1,] -0.95  0.0519      -0.0053       0.1092       0.9025       0.7630
##  [2,] -0.90  0.0338      -0.0034       0.0711       0.8100       0.6848
##  [3,] -0.85  0.0255      -0.0026       0.0535       0.7225       0.6108
##  [4,] -0.80  0.0203      -0.0021       0.0426       0.6400       0.5410
##  [5,] -0.75  0.0166      -0.0017       0.0348       0.5625       0.4755
##  [6,] -0.70  0.0137      -0.0014       0.0288       0.4900       0.4142
##  [7,] -0.65  0.0114      -0.0012       0.0240       0.4225       0.3572
##  [8,] -0.60  0.0095      -0.0010       0.0199       0.3600       0.3043
##  [9,] -0.55  0.0078      -0.0008       0.0163       0.3025       0.2557
## [10,] -0.50  0.0063      -0.0006       0.0132       0.2500       0.2113
## [11,] -0.45  0.0049      -0.0005       0.0103       0.2025       0.1712
## [12,] -0.40  0.0036      -0.0004       0.0077       0.1600       0.1353
## [13,] -0.35  0.0025      -0.0003       0.0052       0.1225       0.1036
## [14,] -0.30  0.0014      -0.0002       0.0029       0.0900       0.0761
## [15,] -0.25  0.0003      -0.0001       0.0008       0.0625       0.0528
## [16,] -0.20 -0.0007      -0.0014       0.0001       0.0400       0.0338
## [17,] -0.15 -0.0016      -0.0034       0.0002       0.0225       0.0190
## [18,] -0.10 -0.0026      -0.0054       0.0003       0.0100       0.0085
## [19,] -0.05 -0.0035      -0.0074       0.0004       0.0025       0.0021
## [20,]  0.00 -0.0044      -0.0093       0.0005       0.0000       0.0000
## [21,]  0.05 -0.0054      -0.0113       0.0005       0.0025       0.0021
## [22,]  0.10 -0.0063      -0.0133       0.0006       0.0100       0.0085
## [23,]  0.15 -0.0073      -0.0153       0.0007       0.0225       0.0190
## [24,]  0.20 -0.0082      -0.0173       0.0008       0.0400       0.0338
## [25,]  0.25 -0.0092      -0.0194       0.0009       0.0625       0.0528
## [26,]  0.30 -0.0103      -0.0216       0.0010       0.0900       0.0761
## [27,]  0.35 -0.0114      -0.0239       0.0012       0.1225       0.1036
## [28,]  0.40 -0.0125      -0.0263       0.0013       0.1600       0.1353
## [29,]  0.45 -0.0138      -0.0290       0.0014       0.2025       0.1712
## [30,]  0.50 -0.0151      -0.0318       0.0015       0.2500       0.2113
## [31,]  0.55 -0.0167      -0.0350       0.0017       0.3025       0.2557
## [32,]  0.60 -0.0183      -0.0386       0.0019       0.3600       0.3043
## [33,]  0.65 -0.0203      -0.0427       0.0021       0.4225       0.3572
## [34,]  0.70 -0.0226      -0.0475       0.0023       0.4900       0.4142
## [35,]  0.75 -0.0255      -0.0535       0.0026       0.5625       0.4755
## [36,]  0.80 -0.0292      -0.0613       0.0030       0.6400       0.5410
## [37,]  0.85 -0.0344      -0.0722       0.0035       0.7225       0.6108
## [38,]  0.90 -0.0427      -0.0898       0.0043       0.8100       0.6848
## [39,]  0.95 -0.0608      -0.1278       0.0062       0.9025       0.7630
## 
## Rho at which ACME = 0: -0.25
## R^2_M*R^2_Y* at which ACME = 0: 0.0625
## R^2_M~R^2_Y~ at which ACME = 0: 0.0528

Here is the summary for imueclt: Country’s cultural life undermined or enriched by immigrants.

summary(sens.cont.4)

## 
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
## 
## Sensitivity Region
## 
##       Rho   ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
## [1,] -0.3 -3e-04      -0.0012        6e-04         0.09       0.0651
## 
## Rho at which ACME = 0: -0.3
## R^2_M*R^2_Y* at which ACME = 0: 0.09
## R^2_M~R^2_Y~ at which ACME = 0: 0.0651

3.2.2 Coefficients of determination

We display the usual graph with one mediator per column, and negative correlations at the top and positive ones at the bottom.

par(mfrow=c(2,4))
plot(sens.cont.1, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens.cont.2, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens.cont.3, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens.cont.4, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens.cont.1, sens.par = "R2", r.type = "residual", sign.prod = "positive", sub="imbleco")
plot(sens.cont.2, sens.par = "R2", r.type = "residual", sign.prod = "positive", sub="imtcjob")
plot(sens.cont.3, sens.par = "R2", r.type = "residual", sign.prod = "positive", sub="imwbcrm")
plot(sens.cont.4, sens.par = "R2", r.type = "residual", sign.prod = "positive", sub="imueclt")

As with the residual coefficient plots, we can see that the effect remains strictly negative for positive products, and that it requires a certain amount of unexplained variance (explained by an unobserved confounder) to get an effect null in reality.

Graphically, one could guess that the corresponding square product for a null effect is around \(R_Y^{2*} R_M^{2*} = 0.1 \times 0.6 = 0.25^2\).

For all mediators except imwbcrm, the ACME is then still guaranteed to be strictly negative as long as the proportion of unexplained variance which can be explained by an unobserved confounder is less than \(25\%\) (for a confounder whose correlations with the mediator and the outcomes are of opposite sign, otherwise it’s okay). We also use the \(\rho\) to state that.

3.3 Conclusions

Our conclusions eventually depend on the mediator:

• imbleco: the proportion of mediated effect is around \(10\%\), the mediation is significant (all the effects are significant), and the result can be considered as robust to the violation of ignorability (see sensitivity analysis),

• imtcjob: the proportion of mediated effect is around \(25-30\%\), the mediation is significant and the result is robust,

• imwbcrm: the proportion of mediated effect is around \(5\%\), with a mediation effect significant at \(90\%\) but not at \(95\%\), and when looking at the sensitivity analysis we can state that the result isn’t robust,

• imueclt: the proportion of mediated effect is around \(25-30\%\), the mediation is significant and the result is robust.

3.4 Discussion

As we can see it by looking at the meaning of those mediations, the main point lies in the causal effect between the treatment and the mediator. One may argue that this relationship hasn’t been studied here, but that would be incorrect as the tests we consider take that relation into account (effect \(a\) with the notations of the report).

Finally, we haven’t included the unused mediators as confounding factors when we studied each mediator in particular (too correlated: that would make the analysis too much tricky), and the unexplained variance in each study could (might) actually be attributed to the unused mediators.

Eventually, even if the proportions of mediated effect can be viewed as low, part of those proportions could actually sum when considering the whole set of mediators.

The two previous points are linked to the issue of multiple mediators raised in the report.

3.5 Bonus: Meaning of the mediations

Let’s have a look at the meanings of the mediations, just to explain the meaning of the results we discuss.

A total effect which is negative (as in our situation) suggests that having minorities living in same current area would make you more inclined to agree with the fact that the Government should be more generous judging applications for refugee status.

The mediation effects we get as results (if proven to be robust enough) suggest that having minorities living in the current area would make people more “immigration friendly” through the idea that:

• immigrants would take less than they put in (imbleco),

• immigrants would create new jobs more than they would take (imtcjob)

• immigrants wouldn’t make country’s crime problems worse (imwbcrm)

• country’s cultural life is enriched by immigrants (imueclt).

In Europe, being in contact with minorities due to a geographical proximity would favor certain positive ideas as listed above, having people more “immigration friendly” (for the results proved as robusts).